The Anisotropic Bak— Sneppen Model 



OO 



o 

B 



-a 
o 

(N 
> 
OO 
(N 
O 
(N 
O 
00 
On 
-i— > 



i 

C 

o 
o 



X 



D A Head* 

Institute of Physical and Environmental Sciences, Brunei University, Uxbridge, Middlesex, UB8 3PH, United Kingdom 

G J Rodgerst 

Department of Mathematics and Statistics, Brunei University, Uxbridge, Middlesex, UB8 3PH, United Kingdom 

(March 9, 1998) 

The Bak-Sneppen model is shown to fall into a different universality class with the introduction 
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I. INTRODUCTION 



The Bak-Sneppen model was originally introduced as 
a crude caricature of biological macroevolution in an at- 
tempt to explain the distribution of extinction sizes ob- 
served in the fossil record Although still widely 
studied in this context, there is also a great deal of inter- 
est in analysing the model from a purely abstract view- 
point. This is because it is currently the simplest and 
most tractable of the class of extremal dynamical mod- 
els, which themselves form a subset of self-organised crit- 
ical systems ||^|||. Extremal dynamical models are so 
called because they are driven by the the selection of 
some globally extremal value which dynamically inter- 
acts with nearby sites. They naturally evolve towards a 
"critical state" (a second order phase transition) without 
any characteristic length or time scales. 

It might appear that the Bak-Sneppen model is well 
suited to adopt the role of the "Ising model" of extremal 
dynamical systems. We believe that this is not the case, 
and in this paper we detail an even simpler version of 
the model which is more open to analysis whilst retain- 
ing all the essential behaviour of the original. The in- 
spiration behind this new model can be most clearly de- 
scribed by analogy with spin systems [||. The Heisen- 
berg spin model is isotropic because the spin vectors 
have no preferred direction. However, when even the 
slightest anisotropy is introduced, a preferred direction 
is created and the system falls into a different univer- 
sality class. Furthermore, this is the same class as the 
highly anisotropic Ising model, where the spin vectors 



can only lie parallel to the direction of quantisation. So 
not only is the Ising model in some sense more general 
than the Heisenberg model, it is also simpler and hence 
more tractable. 

The original incarnation of the Bak-Sneppen model is 
like the Heisenberg model in that it too is isotropic. If 
the analogy with spin systems is to hold true, then the 
introduction of anisotropy into the Bak-Sneppen model 
should result in a different universality class. We find 
that this is indeed the case, at least for one dimensional 
systems, and conclude that, unless there is some rea- 
son for assuming perfect isotropy, it is the anisotropic 
model that should be treated as the general case and 
the isotropic version as a special limiting instance. It 
is also possible to identify a maximally anisotropic Bak- 
Sneppen model which may serve as the true analogue 
of the Ising model for extremal dynamical systems. We 
postpone until Secjv] the question of whether isotropy 
should be assumed in any known application of the 
model. 

This paper is organised as follows. Numerical simula- 
tions of anisotropic systems are described in Sec. [fl] and 
the exponents for the new universality class are given. By 
switching to the multi-trait model, a full solution of the 
maximally anisotropic system is found which explicitly 
demonstrates the crossover to the new class. This solu- 



tion is derived in Sec. Ill alongside the known result for 
the isotropic model. An exact solution for intermediate 
anisotropies was not forthcoming, but by employing an 
alternative means of analysis it is possible to show that 
this new class also applies to any non-zero anisotropy. 
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This is presented in Sec. [V. Finally, in Sec. |y| we dis- 
cuss the applicability of this new class in real situations, 
and consider the potential of the maximally anisotropic 
model in future analytical treatments. 



II. THE ANISOTROPIC MODEL 

Before coming to consider anisotropy we briefly sum- 
marise the isotropic model and some of its known re- 
sults [0J|] . N scalars fi , where i = 1 . . . N , are placed on 
a one-dimensional lattice with periodic boundary condi- 
tions. The fi , known as "barriers," are random numbers 
uniformly distributed on [0,1], although the system be- 
haves in essentially the same manner regardless of the 
particular choice of distribution. At each time step the 
global minimum of all the fi is found, and it and its two 
nearest neighbours are given new random values from the 
same distribution as before. This process is then repeated 
ad infinitum. 

Despite such minimalist dynamics the model exhibits 
a rich variety of non-trivial behaviour. It evolves to- 
wards a statistical steady state in which the bulk of the 
fi are uniformly distributed on [f c , 1], where the thresh- 
old value f c is a function of the lattice dimension and 
connectivity. For the one-dimensional lattice considered 
here, f c ss 0.667. A finite number of barriers form a tail 
on [0, f c ] and it is in this tail that the global minimum is 
always to be found. Both the spatial and temporal cor- 
relation functions are power law in form, signifying the 
existence of a critical state with no characteristic length 
or time scales. The distribution for the absolute distance 
between successive minima Ax is 



behaves in qualitatively the same manner as the isotropic 
model. However, the correlation distributions /jump 



Pjump(Ax) ~ (Ax) 



(1) 



where ir = 3.23 ± 0.02. The probability that the mini- 
mum is at the same site at times to and t + to is given 
by 



-Pall (t) ~ t~ 



(2) 



with tall = 0.42±0.02. This holds true as long as t <C to 
and ageing effects can be ignored [0J|]. 

The model defined above is isotropic because the inter- 
action between the global minimum and the other barri- 
ers is the same in both directions. In other words, if the 
current minimum is fi then barriers /j-i, fi and fi+\ 
are reset, so the minimum is just as likely to jump to 
the left as it is to the right. Consider what happens 
when the rules are altered so that fi—i, fi and fi+2 are 
reset instead. The system now has an inherent bias to 
the right and we would expect an avalanche to be more 
likely to propagate in that direction. This constitutes 
an anisotropic model since there now exists a preferred 
direction for the global minimum to drift. 

We have performed extensive numerical simulations of 
the anisotropic model and have observed that the system 



and Pall have different exponents, 

ALL 



rT 



2.42 ± 0.05 and 



0.59 ± 0.03, so the system falls into a different 
universality class to the isotropic case. Plots of Pall for 
both classes are given in Fig. ^ for direct comparison. 
Pjump is uniformly lower for jumps against the direction 
of anisotropy as for jumps with it, but the same expo- 
nent applies in both directions. The threshold value f c 
also drops, but this is simply due to the increased spread- 
ing out of the avalanche and has nothing to do with the 
loss of isotropy. 

The new universality class is not just restricted to this 
one example. Simulations have shown that if barriers 
fi-a, fi and fi+b are reset, where a and b are arbitrary 
positive integers, then the same class holds for any a =/= b. 
The dynamics can be further generalised by considering 
ranges of sites on each side of the minimum, and either 
selecting all of these sites or just a random sample. Here, 
anisotropy corresponds to a larger range on one side than 
on the other. In all cases the same exponents are found 
for any non-zero anisotropy, although convergence can 
be very slow when the anisotropy is weak, a point that 
will be explicitly demonstrated for the multi-trait model 
in Sec. tVL 
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Logarithm of the number of iterations, ln(t) 

FIG. 1. A log-log plot of Pall(£), the probability for the 
active site to return to its original position after a time t. 
The upper line is from the standard isotropic model and the 
lower line is from an anisotropic system in which the barriers 
fi-i, fi and fi+3 are reset at every time step, where fi is 
the current global minimum. The dashed lines have slopes of 
—0.42 and —0.59, respectively. The data for the anisotropic 
system has been moved upwards to allow for direct compari- 
son with the isotropic case. The simulations were performed 
on an N = 10 4 lattice, for 5 x 10 3 iV iterations in the isotropic 
case and 5 x 10 4 iV in the anisotropic case. 



2 



III. THE MULTI-TRAIT MODEL 

The consequences of introducing anisotropy into the 
Bak-Sneppen model can be more fully investigated by 
switching to the multi-trait framework |j] . In the multi- 
trait model each site has M internal degrees of freedom, 
that is M different barriers rather than just the usual 1. 
At each time step the smallest of all the N x M barriers 
in the system is found and reset. One of the M barriers 
from each of its neighbouring sites is selected at random 
and also reset, so three barriers are reset in total. Then 
the new global minimum is found and the process is iter- 
ated indefinitely. For finite M the system belongs to the 
same universality class as the standard M = 1 model, 
but for M — > oo it falls into a different class and, fur- 
thermore, can be solved exactly. To see why this is so, 
we must first define what is meant by a X-avalanche. 

For any given value of A < f c the global minimum can 
be either greater or less than A. Hence a time series of 
the minimum will consist of regions where it is less than A 
alternating with regions where it is greater than A. Each 
block for which the global minimum is less that A is de- 
fined as a A-avalanche. During a A-avalanche any barrier 
smaller than A is called active since the avalanche can- 
not finish until all of the active barriers have been made 
inactive, that is when they have all been reset to val- 
ues greater than A. There are only two ways in which a 
barrier can get reset, it either becomes the global mini- 
mum or belongs to an adjacent site to the minimum and 
is selected with probability 1/M. However, the latter 
possibility cannot occur in the M — > oo limit since there 
are only a finite number of active barriers in the system 
at any one time, so the probability of selecting one at 
random is vanishingly small. Hence each active barrier 
must eventually become the global minimum and it will 
then initiate a sub-avalanche that can change inactive 
barriers to active, but never the other way around. Fur- 
thermore, since the sub-avalanches from different active 
barriers propagate independently of each other, the ac- 
tive barriers can be reset in any order and there is no 
longer any need to keep track of which is actually the 
global minimum. 

The temporal correlations for M — > oo are the same as 
for a mean field model in which the neighbours of the 
minimum are chosen at random, so the introduction of 
anisotropy will make no difference. Rather than repeat 
the analysis here, we simply quote the main result and 
refer the reader to Q for details of the derivation. If 
P\(t) is the probability that a A-avalanche lasts for ex- 
actly time t, then 

P x {t)~t-*' 2 G(t(\-l) 2 ) (3) 

as A — ► f c = 5, where G(x) is a scaling function that 
tends to a constant value for x — > 0. As expected, Pa 
has the usual mean field exponent of §. Since (||) holds 
independently of the spatial structure of the system, wc 



can already conclude that the threshold barrier value f c 
will be i regardless of the degree of anisotropy. 

Anisotropy will clearly effect the spatial correlations 
and so we present the following analysis in some detail, 
starting with the isotropic model. As in the derivation 
of (||) the algebra is simplified by stipulating that the ac- 
tive barrier always takes the value of 1 when it is reset. 
This makes no qualitative difference to the results. Hence 
the central barrier is always made inactive, but the barri- 
ers reset in each of the adjacent sites may become active 
with probability A. Let g r denote the probability that re- 
setting an active barrier at the origin causes at least one 
of the barriers in site r to become active. Then 1 — g r 
is the probability that no barriers become active, which 
can be related to 1 — g r -i and 1 — g r +\ by the difference 
equation 

1 - g r = (1 - A) 2 + A 2 (l - 5r _!)(l - g r+1 ) 

+ A(l-A){(l- 5r _i) + (l- 9r+ i)} . (4) 

This can be derived by considering what happens when 
an active barrier at the origin is reset. The probability 
of creating no new active barriers is (1 — A) 2 , in which 
case the avalanche will end and site r will definitely not 
become active. This is catered for by the first term on 
the right hand side of (||). Similarly, the second and third 
terms account for the creation of active barriers in one or 
both of the adjacent sites, which may subsequently prop- 
agate to site r with probabilities <? r -i and g r +i, assuming 
g r to be translationally invariant. (^) can be rearranged 
to give 

g r = A(.g,-i + g r +x) - A 2 sv-iflv+i ■ (5) 

If the whole A-avalanche starts from a single active bar- 
rier at r = 0, then g = 1 and (||) can be solved to give 

12 

•9'- = 7 — : 1 ^ ( 6 

(r + 3)(r + 4) 

for A = i , explicitly demonstrating the asymptotic 
power law behaviour g r ~ 1 /r 2 . 

There are many ways in which anisotropy could be 
incorporated into this framework, but for clarity we re- 
strict our attention to just a single definition. At every 
time step the global minimum barrier is found, say in 
site i, and reset. The anisotropic interaction consists of 
randomly selecting one of the M barriers in each of the 
sites i — a and i + b and resetting them both, where the 
parameters a and b are positive integers. Some examples 
are given in Fig. ^. Note that if a and b share a common 
factor, say c, then the system will trivially decouple into 
c independent sublattices. For instance, if a = b = 2 then 
all the even numbered sites will decouple from all the odd 
numbered sites and the two sublattices will evolve inde- 
pendently of each other. Thus we can safely assume that 
a and b are coprime. As a corollary any system with a = b 
is equivalent to the standard model a = b = 1 . Similarly, 
if a is equal to zero we can take 6=1 without loss of 
generality, and vice versa if b = 0. 
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The maximally anisotropic system with a — and 
b = 1 can be solved in much the same way as the isotropic 
case. Since only two barriers are reset at every time step 
anyway there is no need to set the central barrier to 1 as 
before. The resulting difference equation is similar to (Q) 
and can be derived in an entirely analogous manner, 



f 



9r = (l-A) 2 + A 2 (l-.g. r )(l-.g r+1 ) 

+ A(l-A){(l-g r ) + (l-5r+i)} 



(7) 



(a) 









mm 









(b) 



which rearranges to 



9t = 



Aflv-i 



1 - A + X 2 g r 



For A = h this admits the exact solution 



_ J 2Tr fOT T - °' aIld 



for r < 0, 



(8) 



(9) 



so now g r ~ 1/r for large r, giving a power law with an 
exponent of 1. 

An exact expression for Aj^ \ can also be found by 
substituting g r — l/z r into (ra). This gives a linear dif- 
ference equation for the z r , 



_ 1-A 
A 



which can be solved to give 




(10) 



(11) 



For A < z r blows up exponentially in r and so g r 
will exponentially decay to zero. If A > ^ then g r will 
exponentially decay to a constant value for large r, cor- 
responding to a non-zero probability of initiating an in- 
finite avalanche. However, this latter case is of academic 
interest only since the underlying simplification of the 
M — ► oo limit rests on there being only a finite number 
of active barriers at any one time, which is no longer true 



when A > i 



is related to 



the expo- 



The exponent for 
nent for Pjump by n = ra + 1, where g r ~ r~ TR and 
-Pjump(Ax) ~ (Ax) _7r . Hence the analysis given above 
demonstrates that tt changes from 3 to 2 with the intro- 
duction of anisotropy. This should be compared to the 
numerical results in Sec. || for M = 1 systems, where n 
went from 3.23 ± 0.02 to 2.42 ± 0.05. In both cases the 
exponent jumps in the same direction and by a roughly 
similar amount. Furthermore, for M —> oo the exponent 
for Pall (t) ~ t~ TALL obeys t A ll = (2t r )~ 1 . Hence t A ll 
increases from j to 5, and again a similar jump was ob- 
served for M = 1, where tall increased from 0.42 ± 0.02 
to 0.59 ± 0.03. Thus the change in behaviour in the 
M — » 00 limit is also representative of the M = 1 case. 



(c) 

FIG. 2. Examples of various definitions of anisotropy. In 
all cases the solid square is the site with the global minimum 
and the shaded squares are the other sites in which a barrier is 
also reset. The shaded sites are a places to the left of the site 
with the minimum, and b places to its right, (a) The standard 
model a = b — 1. (b) The maximally anisotropic model a = 0, 
6=1. (c) An intermediate case a — 2, b = 3. 



IV. ARBITRARY ANISOTROPY 

It remains to be seen whether systems with interme- 
diate anisotropies do indeed fall into the same univer- 
sality class as the maximally anisotropic model, as im- 
plied by the analogy with the Heisenberg and Ising spin 
models. Unfortunately, the style of analysis adopted in 
the previous section is of little use here since the dif- 
ference equation (|5|) admits no straightforward solutions 
for arbitrary a and b. The difficulty stems from the fact 
that the interactions are now between non-adjacent sites. 
One way around this problem is to find a separate lat- 
tice representation for the avalanche in which only near- 
est neighbours interact. This could then be mapped onto 
the one-dimensional substrate in such a way that nearest 
neighbours on the avalanche lattice map onto interacting 
sites on the substrate. 

To do this unambiguously, it is necessary to employ a 
two-dimensional lattice (n, m) which represents the en- 
tire spatiotemporal extent of the avalanche. The map- 
ping from sites (n, ra) on the avalanche lattice to sites r 
on the one-dimensional substrate is derived as follows. 
The origin (0, 0) corresponds to r = 0. Any given site 
(n, ra) can be reached by taking a steps to the left and 
ra steps to the right, in any order. For arbitrary a and 6, 
the resulting value of r is 

r = aa — rab . (12) 

Each value of r corresponds to the set of points (n^, mi) 
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that obey ([121). Successive points are separated by the 
constant displacement vector 



A„ m = (n i+1 ,m i+1 ) - (ni,m,i) 
= (n i+ i - ni,m i+ i - m l ) 
= (b,a). 



(13) 



To quantify this relationship further, let p nm be the 
probability that site (n, to) is active. Site r will remain 
inactive throughout the entire avalanche only if all of its 
corresponding (n, to) are also inactive, so 



9, 



= 1-11(1 



(14) 



That A„ m is the smallest displacement vector follows 
from the coprime nature of a and b. The mapping 
from (n, m) to r can thus be regarded as a projection 
from the two-dimensional avalanche lattice to the one- 
dimensional substrate. An example is given in Fig. |^(a) 
for the isotropic model. When a ^ b the (n, m)-latticc 
becomes rotated relative to the projection lines, so for 
instance when a — and b — 1 the lattice lies completely 
on its side, as in Fig. ^(b). An intermediate case is given 
in Fig. |(c). 




(a) 

n 



m 



-4-3-2-1 1 2 3 4 



(b) 




(c) 

FIG. 3. The projection from the (n, m) lattice to sites r on 
the one-dimensional substrate. Adjacent sites on the (n, m) 
lattice correspond to sites r that interact, (a) The isotropic 
model a = b = 1. The dotted lines connect all the (n,m) 
that correspond to the same value of r. (b) The maximally 
anisotropic case a — 0, b = 1. (c) An intermediate case a — 2, 
6=1. 



where the product is taken over all the n and to that 
obey ( |l2|) . The advantage of this approach is that varying 
the anisotropy only effects which p nm contribute to (|lj) , 
the p nm themselves are entirely unaltered. Thus a unique 
solution to the p nm exists which, if found, could be ap- 
plied to any anisotropy through (|lj) without modifica- 
tion. 

The next step is to find the solution for the p nm . Each 
site (n, m) has M barriers which, if active, may create 
active barriers in either or both of sites (n + 1, to) and 
(n, to + 1). Since we are still in the M — > oo limit, the 
sub-avalanches initiated by different active barriers are 
independent and can be arranged so as to form a compact 
avalanche on the (n, to) lattice. A site (n, to) can only 
become active if one of its barriers is reset to a value less 
that A due to the interaction with an active neighbouring 
site. The neighbours in question are the two diagonally 
lower sites (n — l,m) and (n,m — 1), so the probabil- 
ity of either event occurring independently is Xp n —i m 
and Ap nm -i, respectively. Thus the difference equation 
for the pnm is 



Pn 



X(Pn- 



Pnr 



Pn — lmPnm- 



1 • 



(15) 



By dropping the second term on the right hand side 
of ( |l5| ) , a linear equation is obtained which has the exact 
solution 



P lin 

III II! 



[n - 



x 



n+m 



(16) 



Except for a missing normalisation factor of 2~( n+m \ 
the coefficients (n + m)!/ (n! to!) describe a binomial dis- 
tribution with equal probability of either outcome. For 
large n and to this binomial is well approximated by a 
Gaussian distribution with mean (n + m)/2 and variance 
(n + m)/4 , 



7r(n + m) 



exp 



1 (n — to) 2 

2 n + m 



(2A) 



(17) 



This expression is vanishingly small except for points 
that lie near the line n — to, which form a non-vanishing 
"tail" . A cut through points of equal n + to shows that 
this tail has a Gaussian cross section of width \\Jn + to , 
so it becomes broader down its length. The behaviour 
down the centre of the tail depends upon the value of A. 
For A i, either blows up or decays exponentially 
according to the factor of (2A) 2ra in (|l7|). At the critical 
point A = 7j this factor becomes unity and instead the 



tail exhibits power law decay p\ 



-1/2 
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Numerical integration of the full difference equa- 
tion (|l5| ) shows that the exact solution of p nm does indeed 
have a Gaussian tail of variance (n + m)/4 which decays 
as a power law for A = i, in agreement with the expres- 
sion for p\" m . However, the exponent for the power law 



decay is different in both cases, p r , 



solution as opposed to p 7 \ 



~ 1 / 2 . The correct expo- 



nent can be recovered by restoring the non-linear term 
in (|l5| ) and instead considering the equivalent continuum 
approximation PQ ]. Let x and y be continuous variables 
corresponding to n and to, and define h(x, y) = p nm - Us- 
ing this notation, 



Vh(x,y) 



dh(x,y) dh(x,y) 
dx dy 

(Pnm - Pn-lm) + {Pnm ~ Pnm-l) , (18) 



and C|l 51) can be rewritten as 



2,2 



-Vh = (2A - l)h - Yh 



(19) 



This can be simplified by making the substitution 
h(x, y) — 1/ z(x, y) and the change of variables u = x + y 
and v — x — y, giving 



(20) 



g = (l-2A)z + A 2 



For A = | the first term on the right hand side of (|2 
vanishes and straightforward integration gives 



z(u, v) 



a 
I 



A(v), 



(21) 



where A(v) is an arbitrary function of v that is found 
from the boundary conditions. There are two sets of 
boundary conditions, one for the line to = and another 
for the line n = 0. It is clear from (|l^) that p n o = A™ 
exactly. The line to = is the same as the line y = 0, 
which maps onto u = v after the change of variables, so 
the first boundary condition is z(u, u) = A~ u . Similarly, 
Pom = A™ and the line n — corresponds to u = —v, so 
the second boundary condition is z(u, —u) = A~". This 
allows for A(v) to be fixed and the full solution is 



z(u, v) 



2 M 



(22) 



Along the tail v — 0, z(u,0) ~ u and so h(x,x) ~ x~~ , 
giving the correct exponent for the decay. However, mov- 
ing away from the tail results in exponential growth in 
z(u, v) and hence exponential decay in h(x, y). Thus the 
continuum approximation predicts the correct exponent 
for the decay of the tail but the wrong cross sectional 
shape, that is, exponential rather than Gaussian. 

The solution for A < | can be found by following ex- 
actly the same procedure. This results in 

((1-2A)A-M + \2) e (i-^)(u-\v\) _ A 2 
z(u,v) = ± — '— . (23) 



As before, this expression increases exponentially away 
from the tail. For v = it reduces to 



z(u,0) 



(1 - A) V 1 - 2 ^" - A 2 
1 - 2A 



(24) 



for the exact which blows up exponentially in u. Hence h(x,x) decays 



to zero when A < \ , giving an exponential cut-off in the 
distribution of avalanche sizes. We note in passing that 
( p3| ) also holds for A > i but, as explained in the Sec. Ill, 
such values of A bear no relevance to actual systems. 

Armed with the solution to the p nm , we can now de- 
rive the exponents of g r for when A = -| . First consider 
the isotropic case a = b = 1. Under the projection in 
Fig. ||(a) all the points down the centre of the tail are 
mapped onto the origin, so g r ^o will only start to receive 
non- vanishing contributions once the tail has become suf- 
ficiently wide. Since the tail broadens like (n + m) 1 / 2 
the first Pnm to contribute to g r will lie on the line 
n + to ~ r 2 , by which point the tail will have already 
decayed to p nm ~ (n + to) -1 ~ r~ 2 . Once these p nm are 
substituted into the infinite product in (p"4|), the leading 
order terms in g r will look something like 



9r 



02 



a 4 



r 6 



(25) 



where the are constants. Hence g r = 0(l/r 2 ), in 
agreement with the exact solution ([|). 

The situation is very different in an anisotropic sys- 
tem a b. As can be seen in Figs. ||(b) and ||(c), the 
tail is no longer vertical but cuts through the projection 
lines at a finite angle, passing over all r > (the preferred 
direction is to the right in both of these examples) . Fur- 
thermore, since the gap between successive p nm mapped 
onto the same r is finite, and the tail broadens with- 
out limit, then for sufficiently large r an arbitrarily large 
number of p nm will contribute to each g r . Each of these 
Pnm will be proportional to 1/r, so the analogous expres- 
sion to ( |25| ) will be g r = 0(l/r) and power law behaviour 
with an exponent of 1 is recovered, in agreement with (|^) 
and numerical simulations. For r < only exponentially 
small pnm contribute and g r <o takes some exponentially 
decaying form. 

Thus the crossover in behaviour from the anisotropic 
to the isotropic model in the M — > oo limit can be 
attributed to the difference between a power law tail 
that moves across the substrate, and one whose cen- 
tre is fixed and can only broaden at a much slower 
rate. The convergence to the new behaviour can be 
very slow, especially for weak anisotropy a w b. Indeed, 
since g r ~ 1/r + 0(1 /r 2 ) the rate of convergence is it- 
self a power law. Slow convergence was also observed in 
the simulations of the M = 1 models is Sec. [n| How- 
ever, for M = 1 the spatial correlations were power law 
in both directions, whereas in the M — > oo limit the cor- 
relations against the direction of anisotropy decay expo- 
nentially. This difference presumably arises because the 
sub-avalanches initiated from different active sites are no 
longer independent for finite M. 
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V. DISCUSSION 



It should come as no surprise that the anisotropic Bak- 
Sneppen model has different critical exponents to its 
isotropic equivalent. Universality classes depend upon 
the dimensionality and symmetries of the model in ques- 
tion, so the loss of symmetrical interactions should result 
in a different class. In spin systems the crossover from 
Heisenberg to Ising behaviour occurs around a given tem- 
perature, which could be very close to the critical tem- 
perature for weak anisotropies. There is no direct ana- 
logue of temperature in the Bak-Sneppen model, where 
the critical state is now the attractor of the dynamics, 
but for weak anisotropies the convergence to the new ex- 
ponents is very slow. Nonetheless, we believe that it is 
the anisotropic class which should now be regarded as 
the general case. 

The isotropic model could still have applications in any 
situation where perfect isotropy can be assumed. The 
question then becomes, do any such situations exist? In 
both of the model's applications we are aware of, we 
think the answer is clearly 'no'. In the biological con- 
text, asymmetry between coevolving species could occur 
for a number of reasons. A graphic example for predator- 
prey relationships is known as the "life/dinner" principle, 
where the asymmetry arises because the prey has more 
to lose from a failed encounter than the predator |fll"| . 
This gets its name from an Aesop's fable, where a dog 
gives up chasing a hare because it is only running for 
its dinner, whereas the hare is running for its life, hence 
"life/dinner" principle [ jl2| . Such asymmetry should re- 
sult in a preferred direction along the food chain, al- 
though the issue is somewhat clouded here by the lack of 
a realistic food web structure (l^ jlj] . A second applica- 
tion of the model has recently been proposed for the pro- 
cess whereby granular materials, such as sand, powder, 
corn flakes etc., settle under perturbations |15|. Here, 
anisotropy would be induced by gravity. 

It was mentioned in the introduction that the max- 
imally anisotropic a — 0, 6=1 system should be more 
open to analysis than the isotropic one. This certainly 
proved to be true in the M — > oo limit studied in Sec. p| , 
where exact solutions were found for all values of A rather 
than just A = \ , as in the isotropic case. The maximally 
anisotropic model may also prove to be more tractable 
in the original M = 1 framework. This claim is not un- 
reasonable and has many precedents. For instance, the 
Zaitsev model is an extremal dynamical system with sim- 
ilar rules to the Bak-Sneppen model, except that a ran- 
dom value is subtracted from the global maximum and 
redistributed equally to its nearest neighbours. Stipu- 
lating that this value is instead only distributed in one 
direction gives rise to an anisotropic variant which can 
be solved exactly, including explicit expressions for the 
critical exponents Another example is provided by 
the abelian sandpile model, where a version in which the 
sand only topples in one direction was solved before exact 



results for the isotropic case were found 17,1a], 

One area that we have not investigated is what hap- 
pens when anisotropy is introduced to lattices with two or 
more dimensions. This opens up the possibility of having 
isotropic interactions parallel to one axis but anisotropic 
interactions parallel to another, the number of permuta- 
tions between the axes increasing with the dimensional- 
ity. Based on the analogy with spin systems, we would 
expect that there would still be just two universality 
classes for each dimension, one for the fully isotropic 
case and one for any non-zero anisotropy. The critical 
exponents for both classes should converge when the up- 
per critical dimension is reached, beyond which the in- 
troduction of anisotropy will make no difference. Work 
is in progress to find at what dimension this conver- 
gence occurs, to help confirm or deny recent claims that 
the upper critical dimension for the Bak-Sneppen model 
is 8 fl|(§. 

Just prior to publication we became aware of a mod- 
ified Bak-Sneppen model by Vendruscolo et al. which 
introduces a preferred direction by a very different mech- 
anism (2l[]. Nonetheless the critical exponents for their 
model appear to match those found for our one, which 
strengthens the case for the universality of the anisotropic 
class. 
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